
#include "vector.h"
#include "matrix.h"

#include <math.h>


/**
 * This assumes the matrix is 2x2. Since this function should
 * not be called unless the data has already been tested by a
 * wrapper no error checking is used. TEST needed...
 */
int _vm_2by2_det(vm_matrix *mtx, double *dest){

  *dest = (mtx->data[0][0] * mtx->data[1][1]) - 
          (mtx->data[1][0] * mtx->data[0][1]);
  return 1;

}

int _vm_3by3_det(vm_matrix *mtx, double *dest){

  vm_vector a;
  vm_vector b;
  vm_vector c;

  a.x = mtx->data[0][0];
  a.y = mtx->data[0][1];
  a.z = mtx->data[0][2];
  a.w = 0;

  b.x = mtx->data[1][0];
  b.y = mtx->data[1][1];
  b.z = mtx->data[1][2];
  b.w = 0;

  c.x = mtx->data[2][0];
  c.y = mtx->data[2][1];
  c.z = mtx->data[2][2];
  c.w = 0;

  vm_vector tmp;
  vm_cross(a, b, &tmp);
  vm_dot(tmp, c, dest);
  //printf("[_vm_3by3_det] %f\n", *dest);
  return 1;

}

int _vm_nbyn_det(vm_matrix *mtx, double *dest){

  /* This recursively takes the determinant of the minors of the */
  /* specified matrix. No error checking is done.                */

  *dest = 0;
  int i;
  for ( i = 0; i < mtx->columns; i++){
    
    double cofact = vm_cofactor(mtx, 0, i);
    cofact *= mtx->data[0][i];
    *dest += cofact;
  }

  return 1;

}

double _vm_rec_nbyn_det(vm_matrix *mtx){

  if (mtx == NULL) return 0;
  if (mtx->rows != mtx->columns) return 0;
  
  if (mtx->rows == 2){
    double dest;
    _vm_2by2_det(mtx, &dest);
    return dest;
  }

  double current_det = 0;
  int i = 0;
  for ( ; i < mtx->columns; i++){
    vm_matrix *minor = vm_minor_matrix(mtx, 0, i);
    current_det += pow(-1, i) * mtx->data[0][i] * _vm_rec_nbyn_det(minor);
    vm_free_matrix(minor);
  }
  //printf("current_det=%f\n", current_det);
  return current_det;

}
